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ABSTRACT 

We present a new method for performing differential emission measure (DEM) inversions on narrow- 
band EUV images from the Atmospheric Imaging Assembly (AIA) onboard the Solar Dynamics Ob¬ 
servatory (SDO). The method yields positive definite DEM solutions by solving a linear program. 

This method has been validated against a diverse set of thermal models of varying complexity and 
realism. These include (1) idealized gaussian DEM distributions, (2) 3D models of NOAA Active 
Region 11158 comprising quasi-steady loop atmospheres in a non-linear force-free field, and (3) ther¬ 
modynamic models from a fully-compressible, 3D MHD simulation of AR corona formation following 
magnetic flux emergence. We then present results from the application of the method to AIA ob¬ 
servations of Active Region 11158, comparing the region’s thermal structure on two successive solar 
rotations. Additionally, we show how the DEM inversion method can be adapted to simultaneously 
invert AIA and XRT data, and how supplementing AIA data with the latter improves the inversion 
result. The speed of the method allows for routine production of DEM maps, thus facilitating science 
studies that require tracking of the thermal structure of the solar corona in time and space. 

Subject headings: Sun: corona - Sun: atmosphere - Sun: activity - plasmas - radiation mechanisms: 
thermal - techniques: observational 


1. INTRODUGTION 


Since the launch of NASA’s Solar Dynamics Observa¬ 
tory (SDO; Pesnell et al. | 2Ql^ in 2010, the Atmospheric 


Imaging As sembly (AIA; Lemen et al. 2012 Boerner 


et al.||2012|) instrument onboard SDO has been deliver¬ 
ing EUV imaging observations of the solar corona with 
an unprecedented combination of temperature coverage, 
spatial resolution, cadence and consistency in data qual¬ 
ity. 

AIA’s simultaneous use of multiple spectral bands 
spanning the range of coronal temperatures also promises 
the ability to diagnose the thermal evolution of the ob¬ 
served systems. However, while the instrument is, by 
design, well-suited for constraining the temperature of 
optically-thin plasma along its line of sight, the routine 
interpretation of AIA data in terms of the temperature 
and density of the emitting plasma remains a difficult 
problem to solve. Because the AIA EUV channels have 
temperature response functions that are generally multi- 
thermal (i.e. they have con tributions from plasma over a 
range of temperatur es, see [Martinez- Sykora et al.||2Qll' 


Boerner et al.||2012|), the thermal distribution of the ob- 
served plasma cannot be directly inferred from the EUV 
images. To learn about coronal temperatures, one must 
separate the multithermal nature of coronal plasma from 
the multithermal response of the EUV channels. 

In the many studies that have used AIA data to at- 
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tempt to infer thermodynamic information about the 
corona, it has often been necessary to use forward fit¬ 
ting or qualitative comparisons with the observed multi- 
spectral characteristics of the plasma. We have devel¬ 
oped a technique to utilize AIA data to their full po¬ 
tential for probing the thermal structure of the corona 
by producing maps of the differential emission measure 
DEM(x, logT, t) at any scale, up to and including the 
full cadence and spatial resolution of the AIA instrument. 
The method we present in this paper has been validated 
against a diverse set of thermal coronal models of varying 
sophistication and realism. 

The remainder of the article is structured as follows. 
In section we formulate the mathematical problem 
underlying DEM analysis and discuss the requirements 
for any inversion method. In section we present our 
novel method for solving the inversion problem. In sec¬ 
tion we test the method against three different classes 
of DEM models. These include log-normal DEMs, distri¬ 
butions from quasi-steady loop models, and distributions 
from a fully-compressible, time-dependent MHD model 
of AR corona formation. In sectionwe apply our DEM 
method to study the thermal distribution of an active 
region observed by AIA. In section we investigate the 
benefits of performing joint DEM inversions using nar¬ 
rowband AIA EUV and broadba nd X-ray imaging data 
from the X-Ra y Telescope (XRTjG olub et al.|[2QC)7|) on¬ 
board Hinode (Kosugi et al.||200 /|). A discussion of the 
implications and scientihc possibilities resulting from this 
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work is given in section 


2. STATEMENT OF THE PROBLEM 
2.1. DEMs from optically thin emission 

Narrowband EUV (or broadband X-ray) observations 
can be related to the physical properties of optically thin 
coronal plasma as an integral over temperature space: 

POO 

Vi = / K,{T) DEM(T)dT, (1) 

Jo 


where yi is the exposure time-normalized pixel value in 
the i-th AIA channel (in units of DN s“^ pixel”and 
Ki(T) is the temperature response function (in units of 


DN cm^ s”^ pixel”^, see Eig.l^.The differential emission 
measure (in units of cm”^ K””) is defined by the relation 
DEM(T)(iT = n‘l{T)dz, where ne{T) is the electron 
number density of plasma at a certain temperature T. 
The integral of DEM(T) over a finite temperature range 
is simply called the emission measure (EM). The aim 
here is to take EUV imaging observations from AIA and 
invert for the emission measure distribution in the so¬ 
lar corona. Although the functions Ki are most sensitive 
to changes in temperature, they also depend on either 
the electron density or pressure. Eurthermore, they de¬ 
pend on the choice of atomic abundances. Th e choices 


we made for this work are given in section 
The DEM inversion problem formulatec. 
example of Eredholm’s integral equation of the first 
kind (see, e.g. Courant & Hilbert|1953[ |Philhps|1962| ) . In 


above is an 


this mathematical framework, the meas urement vector y 
i s an integral transform of DEM(T). [Craig fc Brown 


( 197611^ identified a number of key concerns for DEM 
analysis that result from well-known mathematical prop¬ 
erties of this type of integral equation: 


• Given y^ there may be no solution for DEM(T); 

• Even if a solution exists, it may not be unique; 

• Even if a solution exists, it may be unstable in the 
sense that small changes in y result in large changes 
in DEM(T). This implies measurement errors in y 
may be amplified in the DEM solution; 

• Even if a solution exists, it may not be positive 
(semi) definite. 


These concerns highlight the potential pitfalls of DEM 
inversion analysis. Eor practical purposes, the first may 
not be the most pressing concern. If no solution exists for 
a measurement vector y^ it probably implies the assumed 
physical model is inappropriate. Eor instance, the emit¬ 
ting plasma may not be optically thin, or it may have an 
atomic abundance that is different to that used for com¬ 
puting the response functions Ki{T). Or the plasma may 
have evolved in between exposures in different channels 
such that no single DEM (log T) satisfies Eq. 0 for all 
channels. 

The remaining concerns motivate the following require¬ 
ments for any DEM inversion scheme: 


1. The scheme needs a deterministic way to pick a 
solution out of a family of possible solutions that 
satisfy Eq. Q, and the chosen solution should be 
representative of the emitting plasma; 

2. The solution returned by the inversion scheme 
should be stable. That is, for a noisy measurement 
^ + e, where y represents the noiseless measure¬ 
ment and e represents random errors, the inversion 
scheme needs to return similar DEM solutions over 
different realizations of e; 

3. The solution returned should be positive semidefi- 
nite, i.e. DEM(T) > 0. 

An additional desired property of inversion schemes is 
computational speed. This is especially true for DEM 
analysis of AIA data. In normal operational mode, AIA 
delivers a complete set of seven EUV images every 12 
seconds. This corresponds to ^ 10^ observation vectors 
(^s) per second. Even when subsampling and/or aver¬ 
aging (either spatially or temporally) is used to reduce 
throughput, it would still be desirable for a DEM inver¬ 
sion code to be able to return at least 10^ — 10^ solutions 
per second. 


2.2. Matrix formulation of the problem 

Different DEM inversion methods reported in the liter¬ 
ature have different ways to satisfy (or not) the require¬ 
ments listed in the previous section. It is instructive to 
formulate the DEM inversion problem in matrix form to 
facilitate further discussion. In practice, the integrals in 
Eq. 0 are always approximated as a sum over discrete 
points. Eq. 0 can then be expressed as a set of linear 
integral equations (e.g. see [Craig & Brown 1976). Af¬ 
ter a quadrature scheme is chosen (i.e. in terms of basis 
functions, see Appendix]^, the set of equations can be 
written in the form: 


y = Dx. 


( 2 ) 


The meanings of the matrix elements Dij and compo¬ 
nents of the solution vector Xj depend on the quadrature 
scheme chosen to approximate the integral. 

One choice of a quadrature scheme is the following. 
Let the logarithmic temperature range be divided into n 
neighboring bins. Then Eq. 0 becomes 


Vi = 


n n. 

JlogTj 


logTj+AlogTj 


Ki (log T)DEM(log T)d log T, 


( 3 ) 

where the j-th. temperature bin has range logT G 
[\ogTj,\ogTj -j- AlogTj). Assume that Ki{\ogT) = Kij 
is piecewise constant in each j-th temperature bin, so 


n 

yi = KijEMj , where 
i=i 


EM,= 


I 


log Tj + A log Tj 
log Tj 


DEM(logT)dlogT. 


( 4 ) 

( 5 ) 


^ [Craig & Brown[([l976[) applied their analysis to DEM inversions 
from A-ray line spectra, but the same conclusions apply to DEMs 
from broadband observations. 


Thus Eq. Q has be transformed into the matrix Eq. 0. 
In this quMrature scheme, Xj = EMj so the components 
of the solution vector are simply values of the emission 
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Fig. 1.— Theoretical response functions K{T) for the 94, 131, 
171, 193, 211 and 335 A EUV channels of SDO/ AIA. The respons e 
functions were computed using CHIANTI 7.1.3 ( |Landi et al.|2013] >. 


measure contained in discrete, non-overlapping tempera¬ 
ture bins. K is an m X n matrix with components Kij^ y 
is an m-tuple corresponding to measurements by the AIA 
EUV channels. In practice, we chose to use a set of ba¬ 
sis functions (including both Dirac-delta and Gaussians 
functions) to describe the DEM solution. As detailed in 
Appendix we use a set of Dirac-delta functions and 
three sets m Gaussians of different widths. In total 84 ba¬ 
sis functions are used for a logT grid with 21 grid points 
spanning \ogT/K = 5.5 to \ogT/K = 5.5 at intervals of 
AlogT/K = 0.1. In this case a slightly modified linear 
system is solved but the discussion below still applies. 

Eig. shows temperature response functions for the 
AIA EUV channels as computed by the Solarsoft rou¬ 
tine aia_get-response (/temp, / evenorm, / chian tif ix) 
using the GHIANTI 7.1.3 package ( Landi| 


et d/j |2013|). Goronal abundances specified m 

(a _ co mpilation of 

Grevesse & 


the sun_coronal_ext. abund file 

abundances ___ 

Sanval|1998 Landi et al. |2QQ2p otTIie GHIANTI package 
were used and the pressure was set at = 10^^ 

K cm“^, where k is the Boltzmann constant. The 
chiantifix keyword applies an empirical correction to 
the t emperature resp o nse fu nction of the 94 A chan¬ 
nel (Boerner et al. 2Q14[ ) to account for missing 
transitions in the GHIANTI databaseH 
Although AIA has seven EUV channels, emission in 
the 304 A channel is often optically thick and is not well- 
modeled by GHIANTI under the optically thin assump¬ 
tion. So for AIA DEM analysis the 94, 131, 171, 193, 
211 and 335 A channels are typically used (i.e. m = 6). 
If the number of temperature bins were n = m = 6, 
then one could in principle try to solve for x by multi¬ 
plying both sides of Eg. Q by K~^. This is the solution 
method examined by Craig & Brown (1976). As already 
summarized in section ]2.1[ they pointed out a number of 
problems with this approach. In particular, any noise in 
y will be amplified in x and there is no guarantee that 
the solution will be positive definite. 


^ Older versions also applied an empirical correction to the 
131 A channel but this is no longer the case since version 6 of 
the AIA calibration. 


Eor n > m (i.e. more than 6 temperature bins), Eq. 0 
is under deter mined. This is a well-kn own problem m 
emission measure inversions. Section 12.31 summarizes 
commonly used techniques to deal with this problem. In 
section we present a new method we developed based 
on the concept of sparsity. 

2.3. Methods based on -'m'^nimization 

Most DEM inversion schemes in the literature are 
based on a reduced-y^ approach. That is, the DEM so¬ 
lution is chosen to be one such that 


X^{x) = Y, 

i=l 


I ^Vi j 


(6) 


is minimized. Here 5yi is the uncertainty for the i-th 
channel. This approach is ideal for over determined sys¬ 
tems (i.e. n < m) where it is known that no single 
model will reproduce all n measurements (linear regres¬ 
sion through three or more non-colinear points is one 
example). However, for underdetermined systems it is 
subject to the perils of overfitting. To mitigate this, ad¬ 
ditional constraints are often added to the definition of 
using the method of Lagrange multipliers. The ob¬ 
jective function then becomes + F{x). The regu¬ 

larization term F{x) is used to penalize certain solutions 
with properties that are undesirable. 

A common choice for regularization is the smoothness 
constraint, which is imposed to reject solution s that ex¬ 
hibit oscillatory behavior (e.g. Phillips 1962| Craig & 
Brown||1986l IMonsignori Eossi V Landini||I9911 IHube 


, _ Landlnl||l99llTOg^ 

fc Judge||1995[ [Judge et al.||1997p ~ Smoothness in the 


optimal solution is sought by imposing a penalty term 
F{x) = A "^{xi-i — 2xi + Xi+i)^, where A is the Lagrange 
multiplier and the summand is a finite difference formula 
for the second derivative. This is sometimes called ‘sec¬ 
ond order’ regularization. 

Another choice is the so-called ‘zeroth order’ regular¬ 
ization. In this scheme the penalty term is A||x|| 2 , where 
ll^lb = • X? is the L2-norm (i.e. the Euclidean 

norm) o f x. Zeroth-order regula rizati on was recently 


used b y Hannah & Kontar (2012) and Plowman et al 
([m^ for DEM inve rsio ns of AiA dat^ As d iscussed 
by I Judge et al. (1997) and Plowman et al. (12013), zeroth- 
order regularization has direct correspondence with the 
singular value decomposition (SVD) approach to solv¬ 
ing the underconstrained problem. Both return solutions 
that h ave the minimum total squared EM. It was pointed 
out by Plowman et al. (2013) that zeroth-order regular¬ 
ization ji’THio'irEmTof smoothness constraint since, for 
the same emission measure, a sharply peaked solution 
will have higher total squared EM than a solution with 
a broader distribution. 

Other regularization procedures are available (e.g. 
based on maximum entrop y reg ul arization, see 
Monsignori-Eossi & Landini| |1992| | Judge et HT 
1997p . While regularization helps to mitigate the 


ill-conditioning of the DEM inversion problem, the 
problem remains that x^-minimization schemes gen¬ 
erally do not guarantee solutions that are positive 
definite. Instead, inversion schemes often resort to 


follow-up steps to mitigate negativity (e.g. 

Hannah & 

Kontar||2012 

Plowman et al.||2013). 
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One way to avoid solutions with negative values is to 
perform parametric inversions, i.e. to enforce positive 
definite functional forms of solutions and then to perform 
X^-minimization to find the optimal set of parameters. 
Commonly used fu nctional forms of the DEM function 


data (e.g. time series, i mages, and tomographic magneti c 


include Ganssians (Aschwanden fc Boerner 2011 Guen- 
etar]|2012a|b| ), power laws (|Jordan||l9 76|) 
tion of the two (Gnennon et al.||2013|). 1 


or a com- 
iDiscretized 


bination ^ „ _ 

splines have also been used (llVLonsignori-hossi fc Landini| 
|1992[ [Parent i et al.| 2QQQ[ , see also the xrt_iterative2 
inversion code by Weber, available in the Hinode/XRT 
package in Solarsoft). 

A more comprehensive exploration of the space of pos¬ 
sible DEM solutions can be underta ken using Markov- 
Chain Monte Carlo methods (MCMC; Kashyap & Drake 
1998). These algorithms begin witn a guess at the 
DEM and iteratively apply randomized adjustments in a 
Markov chain, producing a family of DEM solutions that 
can be thought of as a representation of the probability 
distribution function of the actual DEM. The MCMC 
method does not impose a pre-determined functional 
form for the DEM (though it applies some locally vari¬ 
able smoothness based on the shape and coverage of the 
temperature responses), and it is one of the few available 
methods that provides estimates of the uncertainties as¬ 
sociated with the DEM. However, the MCMC method is 
computationally demanding, typically requiring seconds 
to minutes for each observation vector, therefore making 
it not particularly suitable for application to large AIA 
d atasets with high sp atial and temporal resolution. 


Testa et al. 


( 2012) used 3D radiat ive MHD Simula- 


tions of Bifrost ^Guaiksen et al. 2011) to test the relia¬ 


bility of DEM diagnostics applying the MCMC method 
to AIA and Hinode/EIS synthetic data. When applying 
the MCMC method to AIA data they find that though 
the general features and spatial distribution patterns of 
DEMs can be reconstructed, there are some limitations: 
(a) the temperature at which the DEMs peaks is sys¬ 
tematically slightly underestimated; (b) while isothermal 
DEMs are reasonably well reconstructed, inversion solu¬ 
tions for synthetic data generated by broad DEMs, es¬ 
pecially those with significant density structuring along 
the I.O.S., were less accurate. 

3. A NEW METHOD BASED ON SPARSITY 

We propose a new inversion method based on the con¬ 
cept of sparsity, which has received a great deal of atten¬ 
tion in recent years by the compressed sensing commu¬ 
nity. Compressed sensing is concerned with the recovery 
of signals where the number of measurements is less than 
(sometimes much less than) the number of components 
in the reconstructed signals. 

In an underdetermined linear system such as given by 
Eq. the family of solutions satisfying the equation 
resides in an affine subspace of . The challenge is to 
select a solution within this subspace that most faithfully 
represents the underlying scenario. In a series o f papers 
on solu tion s to underde termined linear systems, |Candes| 
& Tao| (e.g. |2006l |2007|) showed that, when compared to 
a least-squares/minimum energy approach, the assump¬ 
tion of sparsity often results in a solution that is a better 
approximation to the real signal. This realization has 
led to immense advances in many fields where the recon¬ 
struction of a linear signal is desired from undersampled 


Donoho|2QQ6 Lnstig et al.j2QQ7 ) 


resonance imaging; see 
Mathematically, the most sparse solution is defined as 
the solution to the optimization problem: 


minimize ||x||o subject to Dx = y. 


(7) 


Here ||:r||o is the LO norm of x, which is just the number 
of non-zero components of x. There is no known efficient 
algorit hm for solving this LQ norm minimization prob¬ 


lem, so Candes & Tao (|200^ instead proposed that one 

dim 


should solve the corresponding LI norm minimization 
problem, namely 


minimize ||x||i subject to Dx = 


( 8 ) 


where ||:r||i = ^ ||^j||- This is the underpinning of our 
i=i 

approach to tackling the EM inversion problem. 

In practice, systematic errors (e.g. in the instrument 
response matrix Kij) and random errors in the measure¬ 
ment vector y means that the sought-after solution may 
not necessarily satisfy Eq. (|^. Eurthermore, for EMs we 
must impose that the solution be positive semidefinite 
(i.e. Xj >0). So our method solves the following linear 
program: 


minimize 

E Xj 

7 = 1 

subject to 

(9) 

Dx 

< 


(10) 

Dx 

> 

max(^— ?/, 0), 

(11) 

X 

> 

0. 

(12) 


The inequality constraint (12) ensures the solutions are 


posi tive semidefinite. The inequality constraints (10) and 
(11) provide some tolerance for the solution to Aviate 
ffoin satisfying Eq. ([^. In this paper, we set the toler¬ 
ance level as r]i = Ci, where is an estimate of the un¬ 
certainty for a pixel count (DN/pixel) in the i-th EUV 
channel divided by the exposure time, can be com¬ 
puted using the function aia_bp_estimate_error, which 
is available as part of the AIA package in Solarsoft. 
aia_bp_estimate_error computes an estimated uncer¬ 
tainty for a given signal level in a given channel based 
on photon counting statistics and a number of instru¬ 
mental effects: read noise, compression and quantization 
round-off, and error in dark subtraction. 

We are unaware of physical principles describing coro¬ 
nal plasmas that would motivate an objective function 
based on the LI norm. However, this choice is appeal¬ 
ing in a number of ways. Eirst of all, this scheme mini¬ 
mizes the number of components (in terms of quadrature 
weights) needed to fit the observations, and in this sense 
it avoids the problem of overfitting. This behavior is con¬ 
sistent with the principle of parsimony (more commonly 
known as Ockham’s Razor). Secondly, this scheme en¬ 
sures positivity of the solution (if a solution is found). 

Thirdly, the problem posed as LPl lends itself to being 
solved by fast numerical techniques. The computational 
requirement of any DEM method is a practical concern 
since AIA delivers data at such a high rate (of order 
10^ observation vectors y per second). The DEM inver¬ 
sion problem posed as LPl is an example of basis pur¬ 
suit. Basis pursuit is a technique commonly employed 
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in the compressed se nsing literature fo r reconstructing 
undersampled signals (jChen et al.|1998|) . Because we re¬ 
quire X > 0, the convex objective tunction ||:r||i reduces 
to the simple linear form Xj. The linear program 
LPl can then be solved effici ently using the simplex al¬ 
gorithm ( Dantzig et aL|p^^ , which is designed to find 
optimal solutions to problems where the objective func¬ 
tion is a linear form and the constraints are posed as lin¬ 
ear inequalities. Our implementation of the DEM inver¬ 
sion code makes use of the simplex function in the IDL 
data analysis package. The implementation of the sim¬ 
plex method in IDL is based on the m ethod as detailed 
section 10. 8 of Nu merical Recipes by [Press, Elanner^ 


& Teukolsky (1986). The computa tiona l speed of the 
inversion code is discussed in section o 
Regardless of the advantages listed above, a DEM in¬ 
version method would be worthless if it only (or mostly) 
returned solutions that are not representative of the 
emitting coronal plasma. In the next section, we present 
results from validation tests of the method. 


4. VALIDATION TESTS 

In this section, we test our inversion method against a 
diverse set of thermal models of varying complexity and 
realism. 


4.1. Gaussian / log-normal DEM distributions 
Log-normal distributions are c ommonly chosen to serve 


as test cases for inversion codes ([Hannah fc KQntar|2Q 12 


Guennou et al. 2012a|b| Plowman et al.||2U13p and as 
functional forms for DEM in versions of AIA data (e.g 


Aschwanden & Boerner 2011). They correspond to Gaus- 
sian functions in log'i' space: 


C(T,Tc,cr) = L^exp 
ayzTT 


(logT-logTe)' 


20-2 


(13) 


where Tc is the peak temperature and a is the Gaussian 
width. The normalization is chosen such that the total 
emission measure is EMq = /o^^dlogTl^ The valida¬ 
tion test for the inversion method was penormed over a 
parameter space of Gaussian distributions ranging from 
log Tc G [5.5, 7.0] and a G [0.0, 0.8]. In the following, the 
total emission measure was set to EMq = 10^^ cm“^. 

Eor each (Tc,cr) model, we computed ^(T) on a tem¬ 
perature grid spanning log T/iG G [5.5, 7.5] at intervals 
of AlogT = 0.0025. Note that for wide DEM distribu¬ 
tions, the total EM contained within this temperature 
range can be significantly (up to 0.2 dex) below 10^^ 
cm“^. We then folded ^ with the AIA temperature re¬ 
sponse functions Ki{T) to generate synthetic count rates 
(DN s“^ pixeD^) for the 94, 131, 171, 193, 211 and 335 
channels. As discussed in section |2.1[ one of the po¬ 
tential pitfalls of DEM inversions is the amplification of 
observational noise in the inversion, which may render 
the solution unrepresentative of the underlying DEM. To 
examine whether this is the case with the sparse inver¬ 
sion method, we generated an ensemble of noisy obser¬ 
vation vectors for each (Tc,cr) model. The magnitude of 


^ Strictly speaking, ^(T) is not the differential emission measure 
as defined in Eq. 0- The two are related by the following relation 
DEM(T) = lnl0T“^^(T). Nevertheless, we will follow the com¬ 
mon practice in the literature and refer to both DEM(T) and C(T) 
as differential emission measure functions. 


the stochastic error e in the AIA channels was estimated 
using the Solarsoft function aia_bp_estimate_error (as 
described in sectionj^for details). For each (Tc, a) model, 
we generated 5000 instances of noisy observation vectors 
Vj + where aj is a random variable drawn from a 

normal distribution. We then performed the DEM inver¬ 
sion on every member of the ensemble. 

For the inversion, we used a temperature grid spanning 
logT/K G [5.5, 7.5], but at a much reduced grid spacing 
of AlogT = 0.1. For a discussion of how the tempera¬ 
ture grid was chosen, see Appendix We have n = 21 
temperature bins for m = 6 AIA channels. We define 
three metrics to quantify the fidelity of the inverted EM 
distribution vs. the input model: 


n 

EM = y]EMj, 

3 


(14) 


log Tem = EM ^ 


EM^- log Tj 

3 


and 


(15) 


w|m = em-i 


y^EMj{\ogTj - logTEM)^ 

3 


.(16) 


They correspond to the zeroth, first and second moments 
of the emission measure distribution. EM^ denotes the 
emission measure contained in the i-th temperature bin 
and EM is the total emission measure, log Tem is EM- 
weighted log temperature and ITem is the effective width 
of the distribution in log T space (indicating the extent of 
multithermality). We use the notation (Q) to denote the 
average of a quantity Q over the ensemble (e.g. (log Tem) 
is the average log Tem taken over the ensemble). 

Figure shows a comparison between the model and 
inverted DEMs. On average, the inversion scheme re¬ 
turns solutions that are representative of the underlying 
input models in terms of the metrics (EM), (log Tem) and 
{Wem)- Over the parameter space of Tc and cr, the code 
is, in general, able to give total EMs with an error of less 
than 10 — 20%, log Tem within an error of 0.21ogT/iG 
and Wem within an error of 0.21ogT/iG. 

Another way to examine the fidelity of inversion solu¬ 
tions with respect to the underlying model is to use joint 
probability density functions (joint PDFs, or 2D his¬ 
tograms). Fig. shows joint PDFs for the three metrics 
EM, Tem and Wem- For all three metrics, the majority of 
the solutions lie along the diagonal, indicating the inver¬ 
sion solutions are generally representative of the under¬ 
lying DEM model. The horizontal error bars indicate the 
uncertainty of a given measurement from the inversion. 
Let P{q^\q^) = P{q^)[f^ P{q^q^)dq]~^ denote the 
conditional probability such that the underlying model 
DEM has metric value q^ given q^ from the inversion. 


Furthermore, let C{q^\q^) = P{q\q^)dq be the asso¬ 
ciated cumulative distribution function. The horizontal 
error bars in Fig. span the range C{q^\q^) = 0.025 to 
C{q^\q^) = 0.975. In other words, given T^j^, one can 
be 95% confident that the underlying model DEM has 
T^ within the range spanned by the er ror bars. 


For t he same Gaussian DEM model, [Guennou et al. 


(2012b) examined the viability of using the six AIA EUV 
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Fig. 2.— Validation test on Gaussian DEM distributions. DEM inversions were carried out using synthetic count rates for six AIA EUV 
channels (94, 131, 171, 193, 211 and 335). Noise was added to the synthetic count rates to generate an ensemble of noisy observation 
vectors for each model. The top, middle and bottom panels show ensemble averages of the three metrics used the quantify the fidelity of 
the DEM solutions (total EM, EM-weighted log T and thermal width) to the underlying DEM model. 
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channels to constrain the DEMs of multithermal plas¬ 
mas. Eor their multithermal DEM inversions, Gaussian 
solutions were assumed with EM, Tc and a as param¬ 
eters. Unregularized x^-minimization (i.e. F{x) = 0) 
was then performed to find the set of parameters that 
best reproduced the synthetic AIA counts. To ensure 
the minimum solution was found for each test case, 
they deployed a GPU-accelerated implementation of the 
inversion code and used a ‘brute force’ approach to 
comprehensively scan the parameter space. Their sys¬ 
tematic study raised concerns about the suitability of 
AIA data for DEM analysis. Eor Gaussian models with 
a = O.llogT/TC, they reported that minimization 
yielded, with high probability, optimal solutions with 
peak temperatures close to the underlying Tc of the 
model. However, the fidelity of the inversions diminished 
with broader DEMs. Eor a = O.TlogT/TC, there was a 
bias toward solutions with Tc = 1 MK for model peak 
temperatures ranging from \ogTc/K = 5 — 7. Also, the 
inverted value for a was only weakly correlated with the 
width of the underlying model for a > 0.3\ogT/K. 

Inspection of Eigs. and [^indicates that the sparse 
inversion method does not suffer from the same s ystem- 
atic effects encountered by Guennou et al. (2012b). Even 
for broad DEMs {a ^ 0.7log'i'/A), the inversion is able 
to return solutions that have, on average, comparable 
width and peak temperatures as the model inputs. Eig.[^ 
shows DEMs inversions for four Gaussian models that are 
warm and hot (log Tc/A = 5.8 and log Tc/A = 6.6, re¬ 
spectively) and narrow and broad (a = O.llogT/A and 
a = 0.71ogT/A, respectively). In each panel, the green 
curve indicates the model Gaussian DEM and the black 
line indicates the ensemble average of the DEM solutions. 
The region bounded by the pair of blue lines represents 
the 95% confidence internal, i.e. 95% of the DEM solu¬ 
tions in the ensemble lie within this region. The ensem¬ 
ble average gives a good representation of the underlying 
input DEM. However, some specific solutions in the en¬ 
semble do contain artifacts. Eor instance, in a Gaussian 
test case with log Tc = 5.8, a = 0.1 (see top panel of 
Eig. 0, some DEM solutions contain a secondary peak 
at \ogF/K > 7.0. 

There is one systematic artifact that is clearly visible 
in the bottom panel of Eig. The underlying model for 
this test case is both hot (logTc/A = 6.7) and broad 
(<j = 0.71ogT/A). At logT/A > 6.6, the ensemble av¬ 
erage has consistently lower emission measure than the 
underlying model. This effect can also be seen in the top 
right panel of Eig. Eor broad DEMs at high temper¬ 
atures, the inversion is consistently underestimating the 
total EM. The middle right panel of Eig.j^shows that, as 
a result, the inversion also systematically underestimates 
Tem when the underlying DEM is hot and broad. As 
the bottom panel of Eig. [^illustrates, the discrepancy is 
mainly due to missing EM in the hig h temperatu r e bins . 
This is consistent with the results of [Testa et al. (2012), 
who also found the temperature systematically iinderes- 
timated for broad DEMs. In section [^ we investigate 
whether using both AIA and XRT observations (the lat¬ 
ter being more sensitive to high temperature plasma) 
helps better constrain the DEM inversions. 


4.2. Nonlinear force-free AR with quasi-steady loops 



5.5 6.0 6.5 7.0 7.5 

Log T/K 



Log T/K 



Log T/K 

Fig. 4. — Comparison between input Gaussian DEMs and in¬ 
verted DEMs from the sparse method for four combinations of Tc 
and (7. In each panel, the green curve indicates the DEM of the 
underlying model. The black solid line indicates the ensemble av¬ 
erage of DEM solutions (for different realizations of noise). The 
pair of blue lines indicates the 95% confidence interval, i.e. 95% of 
the DEM solutions in the ensemble he within the area bounded by 
the blue lines. 


In the second validation exercise, we use DEM distri¬ 
butions from 3D thermal models of NO A A AR 11158 by 
Malanushenko & Schrijver (2015, in preparation). Con¬ 
struction of 3D thermal models of this AR begins with a 
non-linear force-free magnetic field constructed to match 
obser ved AIA loop features (from Malanushenko et al.[ 


2014). The space-filling force-free held is then decom- 
posed into a large number of thin flux tubes (over 7000). 
Elux tubes are defined as the volumes enclosed between 
adjacent field lines, which are traced from a regular grid 








































of seed positions. 

The emission of each flux tube is comput ed individu¬ 
ally assuming a ID qua si-steady atmosphere (Schrijver & 
van Ballegooijeu|2QQ5 ). In each model, the same heating 
scenario is used tor the loop atmospheres of all flux tubes. 
The rendering method used is unique in that it does not 
assume circular cross-sections of flux tubes. Instead, the 
method takes into account the distortio ns experienced 
by the flux tubes along their arc lengths (Malanushenko 
& Schrijver 2013). Given the set of loops for a model, 
the loop temperatures and densities were then resam¬ 
pled onto a Cartesian grid with 2.4 arcsec pixel sepa¬ 
ration (corresponding to 4 pixels for AIA at disk cen¬ 
ter). Integration along vertical columns was performed 
to compute DEM distributions for all pixels. The result¬ 
ing DEM data cube was then folded with AIA response 
functions to obtain the synthetic AIA images shown in 
Fig-IU Although the pixel separation of the synthetic 
images corresponds to 4 AIA pixels, no pixel summing 
has been performed for the synthetic count rates (so this 
is akin to sampling AIA full resolution images at every 
4-th pixel in both spatial dimensions). 

The heating scenarios for the data used here are set 
as follows. The volumetric heating rate e(s) along each 
flux tube is e{s) oc where is the mag¬ 

netic enclosed within a tube, L is the length of the flux 
tube, s is the arc length coordinate and B{s) is the local 
magnetic field strength. The proportionality constant is 
given in terms of the total heating energy input in the 
AR, ^Eh = 5 X 10^^ erg s“^, where the sum is taken 
over all flux tubes. The energy flux entering each flux 
tube is a function of base field strengths, given by 


Eh oc [Bf-V(Bi) + Bt^f{B2) 


(17) 


Following [Schrijver et al. ( [2004 ), we set /(.Bbase) = 
exp{-(Bb ase /500 G)^|. 

We consider two models for our validation exercise. 
Model A is computed with ^ and model B is com¬ 
puted with f3 = 2. Eig. 5 shows synthetic AIA images 
for the two thermal models. These were used as input 
for the sparse inversion code using the same choice of a 
log T grid as for validation test on the log-normal distri- 
bntions . Th e settings for the inversion use d here and in 
section 14.31 are the same as for section 14.11 However we 
do n ot p erform ensemble studies of each pixel as in sec¬ 
tion [T^ and we do not add noise to the synthetic data. 
Rather the purpose is to see how the inversion method 
performs against a diverse set of DEM shapes resulting 
from different physical assumptions (in contrast to the 
idealized log-normal forms considered in the previous sec¬ 
tion). EiguresJ^ and show the comparison between 
the model and inversion for models A and B, respec¬ 
tively. Regions where the comparison is not performed 
are indicated in grey. These are places where the total 
EM is very low (< 10^^ cm“^, too low for inversions given 
AIA’s sensivitity; entire grey region outside of the AR) or 
where the the inversion finds no solutions (isolated pixels 
within the AR). Eor both cases, the inversion is able to 
retrieve values for EM and log Tem that match well with 
the underlying model. As was the case for the valida¬ 
tion exercise using log-normal distributions, the thermal 
width Wem computed from the inversion is less reliable. 
Nevertheless, what is clear from this comparison is that 


the sparse inversion code is able to distinguish the differ¬ 
ent thermal distributions underlying the two models. 

4.3. MHD simulation of Corona Formation in an 
Emerging Active Region 

In the third validation exercise, we use DEM distri¬ 
butions from a time-dependent MHD simulation of the 
formation of AR corona. The simulat i on set up is simi¬ 
lar to the one described in |Chen et al.| (|2014|) except the 
numerical grid spacing is finer. 4'he fully-compressible 
MHD simulation was performed using the Pencil code 


([Brandenburg & Dobler|2002||Bingert & Peter|2011) 
includes a realistic treatment of magnetic held aln 


and 

magnetic held aligned 
thermal conduction in the solar corona. The cartesian 
domain of the simulation spans 147.5 x 73.7 Mm^ in the 
horizontal directions. The bottom boundary is located at 
z = Mm (base of the photosphere) and the top bound¬ 
ary is at 2 : = 50 Mm. At the top boundary the mag¬ 
netic field is matched to a potential field and the mass 
and thermal conductive fluxes are set to zero. Periodic 
boundary conditions were imposed for the side bound¬ 
aries. 

The bottom boundary of the simulation was driven by 
imposing MHD quantities s ampled from an AR forma¬ 
tion simulation described in Rempel & Cheung (2014). 
The latter captures the passage of an untwisted toroidal 
flux rope through the upper 15.5 Mm of the convection 
zone and its eventual e mergence into the overlyin g pho¬ 
tosphere. However the Rempel & Cheung (2014) model 
has a computation dornain that stops at a few hundred 
km above the photospheric base and so does not cap¬ 
ture the corona. By using MHD quantities at the photo¬ 
sphere to set the botton i boundary condit ion for coronal 
simulation, the model of Chen et al. (2014) demonstrated 
how hot coronal loops (at a few million K) spontaneously 
form following enhanced Poynting flux injection at their 
photospheric footpoints. 

Eor this validation exercise, DEMs were computed by 
sampling the simulation cubes of temperature and den¬ 
sity (in a Cartesian domain). Line-of-sight integration of 
the DEM cubes was then performed for a top-down and 
a side view. In both cases, the spatial extent of each inte¬ 
gration column is 432 x 432 km^. This is approximately 
equal to the size of an AIA pixel at disk center. The 
DEM maps were then used for forward synthesis of the 
six AIA EUV images (94, 131, 171, 193, 211 and 335), 
which where t hen used as input for the sparse inversion. 
As in section |4.2[ the primary purpose of the exercise 
here is to test the inversion method against input DEMs 
of different shapes resulting from different physical mod¬ 
els. Here we do not test the inversions on ense mbles 
of noisy synthetic data as we did in section |4.1| Eig- 
ure|^ shows a comparison of total EM, logT and ITem 
between the simulation and inversion results for a top- 
down view. Eig. shows the corresponding comparison 
for a side view. Pixels with total EM < 10^^ cm“^ (in 
the simulation) are omitted from the comparison. The 
comparisons indicate that the sparse DEM method is 
able to retrieve meaningful information about the ther¬ 
mal structure of the active region. In both the top-down 
and side views, the dominant feature is a set of closed 
loops at log Tem^ 6.3. This set of loops connect the 
inner edges of opposite polarity spots in the model AR. 
Loops that emanate from the outer edges of the spots 
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Fig. 5.— Synthetic AIA images (log-scaled) for the two thermal models (top and bottom rows) of NO A A AR 11158. 


are less inclined from the vertical, fan away from the in¬ 
terior of the AR and have XogT^yi/K ^ 5.9. This is 
near the peak of the temperature response function of 
the 171 A channel in AIA (see Fig. [^. The presence of 
these cool peripheral loops in the model (and in the in¬ 
version solution) is consistent with actual 171 A images 
of active regions, which often show bright fans anchored 
at the edge of sunspots fanning away from the interior of 
the AR. An example of such types of structures can be 
seen in Fig. where we applied the inversion code to 
actual AIA ol^ervations of an AR. 

Figure pR| shows four examples of DEM profiles sam¬ 
pled from the MHD simulation. The sampling positions 
are shown as asterisks in Fig. These profiles were 
selected because they have different shapes represent¬ 
ing a diverse range of plasma conditions in the MHD 
model. For instance, the relatively narrow DEM profile 
in Fig. 10 a) is sampled from a cool fan loop with a peak 

b) 

a 


at logT/A 5.75. In contrast, the DEM in Fig 
is sampled near the footpoint of a core AR loop 
peak at \ogT/K = 6.45 accompanied by a broad tail in 
the DEM distribution at lower temperatures. The profile 


shown in Fig. [1Q| ^ 
followed by a oro 


c) also has a peak at \ogT/K = 6.45 
Toad tail. Finally, the profile shown in 
Fig, p^d) has a single peak at \ogT/K = 6.25. 

Recall that the sparse inversion method solves the 


sensitive to the noise level. This implies the details of 
the DEM solution may not be very well constrained (es¬ 
pecially in noisy data). However, it is also clear that the 
inversion can clearly distinguish different types of DEM 
profiles. Whether a DEM inversion (and its solutions) is 
sufficiently good depends on the science question. If the 
aim were to accurately measure slopes of DEM curves in 
the high temperature range (e.g. for constraining coro¬ 
nal heating models), i nversions with AIA data may not 
be the right approach (Warren et al.||2012, also find that 
slopes of DEMs in the high temperature bins have large 
uncertainties). However, as demonstrated in this paper, 
DEM inversions using AIA data can be sufficient for dis¬ 
tinguishing the gross thermal properties between differ¬ 
ent types of coronal environments (e.g. fan loops vs. AR 
core loops). 

5. APPLICATION TO AIA DATA 

The validation exercises in the previous section suggest 
the sparse DEM inversion method is a potentially power¬ 
ful tool for helping us understand the thermal structure 
and evolution of coronal plasma. In this section we pro¬ 
ceed to apply the method to actual AIA observations to 
illustrate its utility. 

Figure shows EM maps from a sparse inversion 
of AIA observations of NO A A AR 11158 and its sur- 


tion maps shown 


effect of pixel averaging (either spatially or temporally) 
on the inversion results can be illustrated by computing 
the effective e^’s consistent with an averaging operation 
(aia_bp_estimate_error supports this) and then using 
those values for the tolerance function of the inversion. 
This effect is illustrated in Fig.[^ which shows inversion 
solutions for three cases: (I) no pixel averaging (dotted 
blue lines), (2) averaging over 9 (red dashed lines) pixels, 
and (3) averaging over 36 (green solid lines). 

The quality of the inversions depends on the signal- 
to-noise ratio of the data under consideration. For this 
reason, one needs to be cautious in order not to over¬ 
interpret inversion solutions. For instance, for each of 


were inverted us- 

tive AR spawned an X2.2 flare (see, e.g. 

Schrijver et al. 

^d assuming no aver- 

2011 

Sun et al. 

2012), two M-class hares and more 


at 2011-02-15122:00. This is more than 20 hours after 
the peak time of the GOES soft X-ray flux during the 
X-flare event. However, AR 11158 continued producing 
a series of C-class flares well after the X-flare. The most 
recent flare that occurred prior to 2011-02-15122:00 
was a C6.6 flare, which began at 2011-02-15119:30 and 
ended at 2011-02-15120:53. 


the cases shown in Fig. 10, the detailed shape (especially 
in the high temperature bins) of the DEM solution is 


The top panel of Fig. 11 shows the total EM contained 
within the temperature range logT G [5.75,6.05]. One 
can clearly identify fan loops anchored at the east and 
west edges of the AR. These cool fan loops are generally 
orien ted away from the AR core (see also [Brooks et al. 
2011). In the temperature range logT G [6.05, 6.65] (the 
next two panels), the high-EM areas mostly delineate 
core loops connecting the leading and following polarities 
of the AR. This transition from cooler peripheral fan 
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Fig. 6. — Validation test on DEM distributions from model A (/3 = 0 in Eq. 0 ) of AR 11158. The layout is the same as Fig.l^ except 
the X- and y-axes here indicate spatial location in the MHD model. Grey pixels indicate regions where the total EM < 10^® cm“^and the 
synthetic counts are too low for DEM inversions (entire area outside of AR) or where the inversion finds no solutions (isolated pixels with 
the AR). The same comparison for a different thermal model (/3 = 2) is shown in Fig.j^ 


loops to warmer core loops is consistent with the t rend 
we saw in the MHD simulation discussed in section lD3l 
In this phase of its life, the AR is still flare productive 
and has high-EM in distinct loops even at temperatures 
above log T/K ^ 6.6. 

When we inspect the DEM maps of the same AR one 
month later (i.e., one solar rotation later), we find very 
distinct differences. At this time, the AR is well into its 
decay phase. Overall, the AR has lower total EM. In 
addition, the EM in temperature bins above log TjK ^ 
6.6 is diminished by at least 2 — 3 orders of magnitude. 
Presumably, the weaker magnetic field strengths and the 
lack of emerging flux and/or shear flows results in less 
energization of the AR complex, leading to lower EMs 
and EM-weighted temperatures. 


The sparse inversion code is fast compared to other 
codes described in the literature for AIA DEMs. The fol¬ 
lowing performance numbers are for a single IDL thread 
running on a MacBook Pro with a 2.6 GHz Intel Core i7 
processor. Eor an inversion problem with six AIA chan¬ 
nels (m = 6), a temperature grid with n = 21 bins and a 
set of / = 4 X 21 = 84 basis functions (see Appendix 
the sparse inversion code computes more than 10^ solu¬ 
tions per second. This does not include the one-time ini¬ 
tialization step of setting up the response matrix, which 
is typically a few seconds and is dominated by disk I/O 
performed by the Solarsoft routine aia_get_response. 
However the start-up time can be amortized as the num¬ 
ber of observation vectors increase. 


5.1. Speed 


6. JOINT AIA & XRT DEM INVERSIONS 
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Fig. 7.— Validation test on DEM distributions from model B (/3 = 2 in Eq. ( |17l ) of AR 11158. The layout is the same as in Fig. 
which shows a similar comparison for a different thermal model (/I = 0) of the AR. 


In our v alida tion tests with log-normal distributions 
(see section 4.1), we found that the sparse method is able 
to return relatively reliable estimates of the total EM, 
EM-weighted logT and (with higher uncertainty) the 
thermal width of the model distribution functions. How¬ 
ever, there was one systematic effect, namely that for hot 
and broad DEM distribution, the method had a tendency 
to underestimate the amount of EM at \ogT/K >6.7 
(see Eig.0. This raises the question of whether the in¬ 
clusion ofan X-ray channel could help improve the DEM 
solution. So we repeated the validation exercise with one 
modification. Instead of just six AIA EUV channels, we 
augment the observation vectors y with an additional 
component corresponding to synthetic count rates ob¬ 
served by the Be-thin channel of XRT. To synthesize 
the XRT count rates we folded the model DEM func¬ 
tions with the temperature response function given by 
the Solarsoft routine make_xrt_temp_resp. This broad¬ 


band X-ray chann el has a response function that peaks 
at \ogT/K = 7.0 dGolub et al.|[2QQ7| ). 

Eignrep!^ shows a comparison of (EM), (logTsM) and 
(ITem) forooth the model and the inversion solutions. 
When we compare this figure with Eig. (showing cor¬ 
responding results for inversions using only six AIA chan¬ 
nels), we find a clear improvement in the fidelity of the 
inversion solutions with respect to the underlying DEM 
distributions. The discrepancy between the inversions 
and the underlying model is much reduced in all three 
metrics. What is interesting is that this improvement is 
not limited to broad and hot DEMs. Even in the parame¬ 
ter regime where the DEMs are centered at lower temper¬ 
atures (say logT < 6.3), the introduction of the Be-thin 
channel results in a marked improvement. Eig. shows 
the corresponding joint PDEs for the joint AlA/XRT 
inversions. Compared with Eig. we see how the intro¬ 
duction of an XRT channel tightens the 95% confidence 




































12 


in 

28.0 

<' 

E 

27.5 

LU 

27.0 

5 

o 

H 

26.5 

26.0 


6.6 


6.3 


6.1 

H 

5.8 

D) 

O 

5.5 


5.3 


5.0 


0.8 

P 

0.6 

D) 

O 

0.4 




0.2 


0.0 


Ground Truth 



X [Mm] 


Inversion 



0 10 20 30 40 50 

X [Mm] 


Inversion - Ground Truth 



0 10 20 30 40 50 

X [Mm] 


0.10 

0.05 


LU 

^ 0.00 




-0.05 I 

-0.10 I 


U) 

o 


0.20 r 

0.10 

0.00 

-0.10 I 

- 0.20 I 


U) 

o 


0.20 r 

0.10 

0.00 

-0.10 I 

- 0.20 I 


Fig. 8.— Validation test on DEM_distri_^tions from an MHD simulation of the formation of coronal loops in an AR (|Chen et al.|2014|). 
The layout is the same as in Figs. and Within the AR core loops there are regions where the relative error in total EM logElvl) 
can be of order unity (black pixels in the top right panel). This is mostly restricted to regions where total E M < 10^^ cm“^. The same 
comparison for a side view is shown in Fig. DEM profiles for positions marked by asterisks are shown in Fig. |10| 
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Fig. 9.— Validation test on DEM distributions from an MHD simulation of the formation of coronal loops in an AR ( |Chen et al. poTil i. 
The same comparison for a top-down view is shown in Fig. 
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Fig. 10. — Examples of DEM profiles sampled from the MHD 
simulation at the four positions marked by asterisks in Fig.[^ The 
black solid lines indicate the underlying DEMs from the simulation 
(i.e. ground truth). In each example, three solutions are shown, 
corresponding to inversions with tolerance levels with no pixel av¬ 
eraging (blue dotted lines), and for averaging over 9 (red dashed 
lines) and 36 pixels (green solid lines). 
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Fig. 11.— DEM maps of NOAA AR 11158 and its surround¬ 
ings obtained by the sparse inversion method. The six AIA EUV 
images were taken at 2011-02-15722:00. The field-of-view spans 
1200x480 arcsec^ and is centered at {x^y) = (600,-268) arcsec. 
The color-coding indicates the total EM contained within a logT 
range indicated in the bottom left corner of each panel. 


intervals for all three metrics. 

This exercise demonstrates two important points. 
Firstly, the sparse DEM inversion method can easily be 
extended to take into account data from instruments 
other than AIA. Secondly, the inclusion of X-ray data 
(such as from XRT) can help improve the DEM inver¬ 
sion results (see also Hanneman & Reeves 2014 for a 
discussion of the impact of performing joint AlA-XRT 
DEM inversions). In future work, we will perform joint 
AIA-XRT inversions on real data. This requires a careful 


examination of the intercalibration between the instru¬ 
ments and is outside the scope of the present study. 

7. DISCUSSION 

By delivering full-disk EUV observations of the Sun 
at high cadence, spatial resolution and regularity, 
SDO/AIA has so far proven immensely valuable for stud¬ 
ies of the dynamical behavior of coronal features in EUV 
images. The interpretation of optically thin EUV fea¬ 
tures in AIA images in the thermal domain requires 
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Log Emission Measure [cm"] 


26.00 26.50 27.00 27.50 28.00 28.50 


EM in log T/K= [5.75,6.05] 


EM in log T/K= [6.05,6.35] 


EM in log T/K= [6.35,6.65] 


EM in log T/K= [6.65,6.95] 


EM in log T/K= [6.95,7.25] 


2011-03-15122:00 


Fig. 12.— Same as Fig. |11| but one month later, at 
2011-03-15722:00. At this time, AK 11158 has decayed substan¬ 
tially and there is very little EM at temperatures above log TjK = 
6.65. 


either DEM analysis and/or sophisticated MHD mod- 
eling with forward synthesis of observables (e.g. |Peter 

et al.|2004l Guc 

iksen & Nordlund|2005l 

Mok et al. 1*20051 

Hansteen et al. 

120101 IPeter & Bingert 

20121 |Jin et al. 

20131 |van der E 

Lolst et al. 20141 Chen et al.||2014p. The 


ing of the physical processes operating in solar plasma. 
Forward modeling of synthetic observables from numer¬ 
ical models and a direct comparison of these computed 
quantities with observation data provides the most direct 
method for testing the validity of the models. However, 


since numerical simulations are computationally expen¬ 
sive, they generally sample a very limited range of param¬ 
eter regimes and physical scenarios. Thermal analysis 
techniques such as DEM inversions are therefore impor¬ 
tant for probing the parameter range available to solar 
phenomena. Using the concept of sparsity, we have im¬ 
plemented a high throughput DEM inversion code that 
returns positive semidefinite solutions with a computa¬ 
tional speed that is suitable for producing DEM maps 
from AIA images. 

Before applying any DEM method to actual solar data, 
one must validate the method against a variety of model 
thermal distributions to ensure the method returns so¬ 
lutions that are representative of the physical scenario 


dation tests (see also iTesta et al. 

2012 Guennou et al.| 

2012a|b Plowman et al.||2013|) on 

three classes ot DEIVI 


4.1), DEMs from quasi-steady loops loa ded o n non-linear 


force-free field models of an AR (section 4.2), and DEMs 
from a fully time-depen dent MHD simulation of AR 
corona formation (section |4.3[ ). By testing the inversion 
method against models of varying complexity based on 
different physical assumptions, we mitigate the tendency 
to tune the inversion method to one particular type of 
DEM distribution (a type of overfitting). 

To illustrate the utility of our new inversion method, 
we applied it on DEM inversions of AIA observations of 
NOAA AR 11158 at two different phases of its life. DEM 
maps of the AR during its flaring phase (Eig. [IT] ) show 
high-EM core loops above 3.5 million K. In comparison, 
DEM maps of the AR during its decay phase show a 
clear deficit of material at similar temperatures. In both 
phases, high-EM structures in the low temperature bins 
(0.5 — 1 million K) outline fan loops predominantly an¬ 
chored at the periphery of the AR. 

In section we showed how supplementing AIA im¬ 
ages with an XRT channel (e.g. Be-thin) improves the 
quality of sparse DEM inversions. This hypothetical ex¬ 
ample serves to illustrate how the sparse DEM method 
can easily be adapted and/or extended for inversion of 
data from other (or multi-) instruments. 

Whether DEMs with AIA data are appropriate for 
probing the thermal structure of the solar corona de- 
pendson the science question. As discussed in sec¬ 
tion |4.3[ the slopes of DEMs in high temperature range 
are not well constrained. So if this measurement were 
crucial for testing coronal heating models, the present 
approach is perhaps not appropriate. However, in this 
paper we have demonstrated the ability of the sparse 
inversion code (as applied to AIA data) to yield DEM 
solutions that allow one to distinguish between different 
types of thermal structures in the solar corona (e.g. cool 
fan loops vs. AR core loops). The present work opens 
up the possibility of routine production of DEM maps 
using AIA data (from a validated method) even in the 
absence of complementary coverage from XRT and EIS. 
This allows AIA to fulfill its promise to help researchers 
probe the thermal distribution and evolution of the solar 
corona. 


The authors acknowledge support from NASA’s 
SDO/AIA contract (NNG04EA00C) to LMSAL. AIA is 
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Fig. 13.— Similar to Fig. [^except here we performed the DEM inversions on observation vectors comprising six AIA EUV channels (94, 
131, 171, 193, 211 and 3355 and the Be-thin channel from XRT. Noise was added to the synthetic count rates to generate an ensemble 
of noisy observation vectors for each model. The top, middle and bottom panels show ensemble averages of the three metrics used the 
quantify the fidelity of the DEM solutions (total EM, EM-weighted log T and thermal width) to the underlying DEM model. 
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Fig. 14.— Similar to Fig. [^except here we performed the DEM inversions on observation vectors comprising six AIA EUV channels (94, 
131, 171, 193, 211 and 335) and the Be-thin channel from XRT. Compared to using AIA data only, the joint AIA-XRT inversion provides 
DEM solutions that better match the underlying model. 
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APPENDIX 


QUADRATURE SCHEME 

Let i = 1, 2,m denote the index over a set of wavelength band channels and/or line spectra. Let the DEM function 
be written in terms of a set of positive semidefinite basis functions {bj{\ogT) >0 | /c = 1, 2,/}, viz. 


i 

DEM(logr) = J2h{logT)xk, (Al) 

k=l 

with quadrature coefficients Xk >0. Approximating the integrals in equation 0 as sums in log T space, we have 

n I 

KijBjkXk^^ogT, (A2) 

i=l k=l 


where j = 1,2, ...,n is the index over temperature bins, Kij = iC^(logTj) and Bj^ = 6/c(logTj). The response matrix 
K = (Kij) has dimensions m x n. The basis matrix B = has dimensions n x I, with the k-th column vector 

corresponding to the k-th basis function bk(iogTj). Defining the dictionary matrix D = KB, the set of integral 
equations 0 can be written in matrix form: 

y = Dt, (A3) 

where the sought-after solution vector x is an /-tuple with components x/^AlogT (k = 1,2, ...,/ ). When the number 
of basis functions exceeds the number of image channels (i.e. I > m), the linear system Eq. ( |A3[ ) is underdeter mined. 

Eor the results in this paper, we use an equidistant grid in logT with AlogT = 0.1 ranging from logT = 5.5 to 7.5 
(i.e. n = 21). Over this temperature grid, the set of Dirac-delta basis functions \ k = ...,n} is 

6D--(log Tj) = 1, if log Tj = log Tk, (A4) 

= 0, otherwise. (A5) 


Recall that the basis matrix B consists of column vectors corresponding to basis functions. So for the set of Dirac-delta 
functions = I (the identity matrix). 

In addition to Dirac-delta functions, we also use basis functions consisting of truncated Gaussians. Each Gaussian 
function of width a generates a set of basis functions {b^ | /c = 1, where 


bki^ogTj)=^^P 


(logT,- - logTfe)2 

n2 


if I logT,- - logTfcl < 1.8a. 


(A6) 


= 0, otherwise. 


(A7) 


The Gaussian basis functions are truncated (i.e. set to zero) for values of log Tj outside the temperature grid used 
for inversions. The corresponding basis matrix for this set is denoted B^. Different sets of basis functions can be 
combined by concatenating their associated basis matrices. Eor the inversions shown here, we use the combined basis 
matrix 


B = (B^ 


Dirac |ga= 0 . 1 1-[50=0.2 |-[50=0.6 


B“ 


|B“ 




(A8) 


Note the individual Gaussian basis functions are not normalized by their sums (i.e. all have maximum value of unity at 
their peaks). So given multiple solutions that equally fit the data, the method will prefer a solution consisting of a single 
broad Gaussian over solutions consisting of multiple narrow Gaussians (and/or Dirac-delta functions). Empirically, 
we find the choice of not normalizing the Gaussian basis functions results in better inversions results (more on this 
below). Because n = 21, B as indicated above (see Eig. 15 for a graphical representation) has dimensions 21 x 84 
and D has dimensions m x 84. With six AIA channels, m = 6. Even when AIA is augmented by XRT or EIS data. 


m 84. This makes Eq. (A3) a highly underdetermined system, which we solve by the method of basis pursuit (see 
section . 

Seeking to solve this underdetermined system is the same as the following geometric problem. Suppose we aim to 
express some given column vector y (of dimension m) as the linear combination of members drawn from a family of 

column vectors. Let this family of column vectors be denoted {dk \ k = 1,...,/}. The goal is to find coefficients Xk 


such that y = X]/c=i ^kdk^ which is equivalent to the linear system (A3) if the dictionary matrix D is constructed by 

concatenating the d/^’s side by side. Because / > m, {d^ | /c = 1,...,/} is an overcomplete set of po ssible basis vectors 
(i.e. dictionary) for building up y. So the non-uniqueness of a DEM solution satisfying Eq. (A3) is the same as the 
multiplicity of ways to find a basis for y. Basis pursuit addresses this by seeking a solution that minimizes the L 1 -norm 
\xU. In oth er words, basis pursuit finds the most sparse representation of y from an overcomplete dictionary D (Ghen 
et al.|p^M ). 
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Fig. 15.— Each column vector in the basis matrix B given in Eq. ( |A8| > corresponds to a basis function. The basis matrix shown in 
the upper panel results from four sets of basis functions. The leftmost set corresponds to Dirac-delta functions and the remaining three 
correspond to truncated Gaussians of width a = 0.1, a = 0.2 and a = 0.6 log T/iF. In the shading used above, black indicates unity and 
white indicates zero. The lower panel shows the likelihood that a basis funct ion (i.e. a column in the basis matrix) has a non-zero coefficient 
in the DEM inversion test performed on log-normal DEMs (see section [4.1| . The black line corresponds to the case where the basis matrix 
shown in the upper panel (i.e. Dirac-Delta functions and three sets of llaussians) is used. The red line corresponds to the case when only 
Dirac-delta functions are used as the basis. 


It is worth comparing the fidelit y of the inversions for different bases. For the results shown throughout this paper, 
we used the basis given by Eq. (A8). This was a choice made after performing validation e xerci ses for a number 
of different bases. For example, consider the validation exercise on log-normal DEMs (section |4.1[) . Fig. ^ shows a 
comparison between the model DEMs and the inversion solutions in t erms of total EM, EM-weighted log T and thermal 
width. If we choose the basis B = = I, the linear system (A3) reduces to Eqs. (|^ and ([F). For this choice of B, 

the least Ll-norm principle has a direct physical counterpart, namely a solution is sought suoi that the total EM is 
minimized. While the connection with a physical principle is appealing, this choice clearly gives an inferior inversion 
resul t (se e Figs. 

Eq. (A8) but wr 

scheme chooses to express the solution as a sum of isothermal components). This motivated us to use Gaussian basis 
functions that have maximum values of unity regardless of their widths. This effectively introduces a preference for 
broad solutions over narrow DEM solutions. 

There are likely better choices of basis functions for DEM inversions and the optimal set of basis functions (if it 
exists) may ultimately depend on the problem of interest. How to choose optimal choices of basis functions is an open 
question but theory and perhaps numerical simulations can provide guidance for future improvements. 


I6|and I7|). The following point is also worth noting: Inversions performed usir^the basis given by 
noruiauzed Gaussians result in the same results as indicated in Figs. 16 and[l7| (i.e. the inversion 


DEPENDENCE OF DEM INVERSIONS RESULTS ON THE CHOICE OF THE RANGE FOR logT 

The DEM test results presented in this paper use a logT grid that spans \ogT/K = 5.5 to \ogT/K = 7.5 at 
intervals of AlogT/iC = O.I. As Fig. [^indicates, however, some EUV channels of AIA (e.g. 335) also show significant 
response to plasma at \ogT/K <5.5. Our decision to restrict the DEM inversion to logT/iU = 5.5 is motivated by 
tests in which we varied the lower bound of the temperature grid. F or ex ample, Fig. p!8] shows the spatially-averaged 
DEM distribution for model A of the test cases considered in section |4.2| (quasi-steady loop models by Malanushenko 
et ah). For generating synthetic AIA images (shown in the top row of Fig. [^, we used the model DEM values 
between \ogT/K = 5.0 to logT/iU = 7.0 (though there is no plasma in the model with logT/iU > 6.8). We then 
performed pixel-by-pixel sparse DEM inversions on the synthetic AIA data and computed the spatial average of the 
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Fig. 16. — Similar to Fig. [^but for a quadrature scheme that uses only Dirac-delta basis functions. The quality of the DEM inversions 
is clearly inferior in this case. The black patches in the top right corners of the plots means the inversion is returning total EMs and 
EM-weighted log T values with discrepancies exceeding the range displayed by the color table. 


DEM solutions. 

The black line in Eig. 18 shows the underlying model DEM. The green line shows the DEM solution for a \ogT/K 
grid spanning 5.5 to 7.5, with AlogT/K = 0.1. The blue and red lines show the corresponding solutions when the 
lower bound of the temperature grid is \ogT/K = 5.2 and \ogT/K = 5.0, respectively. Between \ogT/K « 5.6 
and \ogT/K ^ 6.3, all three solutions closely track the underlying model DEM. However, the blue and red solutions 
clearly show spurious enhancements of EM (relatively to the underlying model) at transition region temperatures 
(logT < 5.5). Associated with this enhancement in the transition region is a deficit of EM between \ogT/K = 6.4 and 
logT/K = 6.6 (corresponding to the cores loops of the model AR). Erom this type of test, we decided to restrict the 
temperature grid to logT/K = 5.5 and above for the validation exercises described in the paper. 
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Fig. 17.— Similar to Fig.j^but for a quadrature scheme that uses only Dirac-delta basis functions. The quality of the DEM inversions is 
clearly inferior in this case. Here the inversion solutions are expressed as the sum of a few isothermal components and result in distributions 
that are much more spiked than the underlying DEM distributions. 
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Fig. 18. — Dependence of DEM solution with respect to choice of temperature grid. Plotted above are spatially-averaged DE M dis¬ 
tributions of the model AR shown in Fig. The black line indicates the underlying DEM model (model A, see section |4.2| ). The 
green line indicates the sparse inversion soluOon for a logT/AT grid spanning 5.5 to 7.5. The blue and red lines indicate solutions where 
the lower limit of the temperature grid is set to logT/AT = 5.2 and logT/AT = 5.0, respectively. Of the three solutions, the one with 
logT/AT = 5.5 as the lower temperature limit (green line) shows the best match against the underlying model DEM over the temperature 
range logT/AT G [5.5, 6.7]. 






